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Abstract. It is well documented that a model for the underlying asset price process that seeks to capture the 
behaviour of the market prices of vanilla options needs to exhibit both diffusion and jump features. In this paper 
we assume that the asset price process 5 is Markov with cadlag paths and propose a scheme for computing the 
- - - law of the realized variance of the log returns accrued while the asset was trading in a prespecified corridor. We 

, thus obtain an algorithm for pricing and hedging volatility derivatives and derivatives on the corridor-realized 

I variance in such a market. The class of models under consideration is large, as it encompasses jump-diffusion 

and Levy processes. We prove the weak convergence of the scheme and describe in detail the implementation of 
the algorithm in the characteristic cases where S is a CEV process (continuous trajectories), a variance gamma 
process (jumps with independent increments) or an infinite activity jump-diffusion (discontinuous trajectories 
with dependent increments). 



o 

(N 

; 1. Introduction 

d 

If^ I Derivative securities on the realized variance of the log returns of an underlying asset price process trade 

O^i actively in the financial markets. Such derivatives play an important role in risk management and are also 

used for expressing a view on the future behaviour of volatility in the underlying market. Since the liquid 

^ . contracts have both linear (variance swaps) and non-linear (square-root = volatility swaps, hockey stick = 

^ ■ variance options) payoffs, it is very important to have a robust algorithm for computing the entire law of 

' the realized variance. Often such contingent claims have an additional feature, which makes them cheaper 

• ■ and hence more attractive to the investor, that stipulates that variance of log returns accrues only when the 
in ■ , , 

O ' spot price is trading in a contract-defined corridor (see Subsection 12.11 for the precise definitions of such 
^ I derivatives). 

It is clear from these definitions that, in order to manage the risks that arise in the context of volatility 
^ I derivatives, one needs to apply the same modelling framework that is being used for pricing and hedging 
vanilla options on the underlying asset. It has therefore been argued that the pricing and hedging of volatil- 
ity derivatives should be done using models with jumps and stochastic volatility (see for example lITOl . 
Chapter 11). In this paper we propose a scheme for computing the distribution of the realized variance and 
the corridor-realized variance when the underlying process S = {S)t>o is a Markov process with possibly 
discontinuous trajectories, thus obtaining an algorithm for pricing and hedging all the payoffs mentioned 
above. Our main assumption is that the Markov dimension of S is equal to one (i.e. we assume that the 
future and past of the process S are independent given its present value). We do not make any additional 
assumptions on the structure of the increments or the distributional properties of the process S. This class of 
processes is large as it encompasses one dimensional jump-diffusions and Levy processes. 



We would like to thank Martijn Pistorius for many useful discussions. 

1 



2 



HARRY LO AND ALEKSANDAR MIJATOVIC 



The algorithm consists of two steps. In the first step the original Markov process S under a risk-neutral 
measure is approximated by a continuous-time finite state Markov chain X = (X,),>o. This is achieved 
by approximating the generator of 5 by a generator matrix for X. The second step consists of pricing the 
corresponding volatility derivative in the approximate model X. It should be stressed that the two steps are 
independent of each other but both clearly contribute to the accuracy of the scheme. In other words the 
second step can be carried out for any approximate generator matrix of the chain X. In specific examples 
in this paper we describe a natural way of defining the approximate generator matrix (see Section |4] for 
diffusions and Section [5] for processes with jumps) which is by no means optimal (see monograph 111 
for weak convergence of such approximations and |[T6l for possible improvements) but already makes the 
proposed scheme accurate enough (see the numerical results in Section [6l). 

In the second step of the algorithm we approximate the dynamics of the corridor-realized variance of the 
logarithm of the chain Z (i.e. the variance that accrued while X was in the prespecified corridor) by a Poisson 
process with an intensity that is a function of the current state of the chain X. This approximation is obtained 
by matching G N instantaneous conditional moments of the corridor-realized variance of the chain X. This 
is a generalisation of the method proposed in HI, which in the framework of this paper corresponds to ^ = 1 
and only works in the case of linear payoffs on the realized variance (i.e. variance swaps) as can be seen 
in Tables [51 [6] and |7] of Section |6] Using k strictly larger than one improves considerably the quality of the 
approximation to the distribution of the corridor-realized variance for S. In fact if 5 is a diffusion process, 
then our algorithm with k = 2 produces prices for the non-linear volatility payoffs (e.g. volatility swaps and 
options on variance) which are within a few basis points of the true price (see Tableland Figure [Tell- If the 
trajectories of S are discontinuous, then the scheme with k = 3 appears to suffice (see Tables |6] and [T] and 
Figures |2b] and [3all. Note also that in [14] we provide a straightforward implementation of our algorithm 
in Matlab for ^ = 3. Furthermore in Section |3] we prove the weak convergence of our scheme as k tends to 
infinity (see Theorem 13. II ). 

The general approach of this paper is to view continuous-time Markov chains as a numerical tool that is 
based on probabilistic principles and can therefore be applied in a very natural way to problems in pricing 
theory. It is worth noting that there is no theoretical obstruction for extending our scheme to the case when S 
is just one component of a two dimensional Markov process (e.g. S is the asset price in a stochastic volatility 
model) by using a Markov chain to approximate this two dimensional process. The reason why throughout 
this paper we assume that S itself is Markov lies in the feasibility of the associated numerical scheme. If 
S is Markov the dimension of the generator of X can be as small as 70, while in the case of the stochastic 
volatility process we would need to find the spectra of matrices of dimension larger than 2000. This is by 
no means impossible but is not the focus of the present paper. 

The Uterature on the pricing and hedging of derivatives on the realized variance is vast. It is generally 
agreed that either the assumption on the independence of increments or the continuity of trajectories of the 
underlying process needs to be relaxed in order to obtain a realistic model for the realized variance. In 
the recent paper Q model independent bounds for options on variance are obtained in a general continu- 
ous semimartingale market. The continuity assumption is relaxed in |7 |, where a class of one dimensional 
Markov processes with independent increments is considered and the law of the realized variance is ob- 
tained. A perfect replication for a corridor variance swap (i.e. the mean of the corridor-realized variance) in 
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the case of a continuous asset price process is given in ||6l. For other contributions to the theory of volatihty 
derivatives see HI and the references therein. The main aim of this paper is to provide a stochastic approxi- 
mation scheme for the pricing and hedging of derivatives on the realized (and corridor-reaUzed) variance in 
models that violate both the above assumptions, thus making it virtually impossible to find the laws of the 
relevant random variables in semi-analytic form. 

The paper is organised as follows. Section|2]defines the approximating Markov chains and gives a general 
description of the pricing algorithm. In Section [3] we state and prove the weak convergence of the proposed 
scheme. Section|5](resp.|4l) describes the implementation of the algorithm in the case where the process S is 
an infinite activity jump-diffusion (resp. has continuous trajectories). Section |6] contains numerical results 
and Section |7] concludes the paper. 

2. The k conditional moments of the realized variance 

Let S = {St)t>o be a strictly positive Markov process with cadlag paths (i.e. each path is right-continuous 
as a function of time and has a left limit at every time t) which serves as a model for the evolution of the risky 
security under a risk-neutral measure. Note that we are also implicitly assuming that 5 is a semimartingale. 

2.1. The contracts. A volatility derivative in this market is any security that pays 0([log(5')]7-) at maturity 
T, where <p : M+ ^ M is a measurable payoff function and [log(5)]r is the is the quadratic variation of the 
process log(5) = (log(5;));>() at maturity T defined by 

/ \ ^ 

(1) [log(5)]r := lim ^ 

/7 — ;.oo 

t"en„,i>i 

where n„ = {t^j",. . . ,t"}, « S N, is a refining sequence of partitions of the interval [0, T]. In other words 
f« = 0, = r, n„ C n„+i for all « G N and hm„^oomax{|ff - J : / = 1, . . . ,«} = 0. It is well-known 
that this sequence converges in probability (see |[T3l . Theorem 4.47). Many such derivative products trade 
actively in financial markets across asset-classes (see lUl and the references therein). 

A corridor variance swap is a derivative security with a linear payoff function that depends on the accrued 
variance of the asset price S while it is trading in an interval [L,U] that is specified in the contract, where 
<L <U <oo. More specifically if we define a process 

(2) S,:=max{L,min{St,U}}, VfG[0,oo), 

then for a given partition n„ = {fg,?", . . . of the time interval [0,7] the corridor-realized variance is 
given by 

2 

/ V „ \ 

(3) l[L,C/](%,) + l[L,C/]('^r,")-l[L,f/](%i)l[L,f/]('^f,") 

f,"en„,r>i 

where 1[l.u] denotes the indicator function of the interval [L,U]. In practice the increments tf — t" ^ ususally 
equal one day. The square bracket in the sum in Q ensures that the accrued variance is not increased when 
the asset price 5 jumps over the interval [L, 

The one sided corridor-realized variance was defined in [4]. Definition (1.1) in [4] (resp. (1.2) in ll4l ) 
corresponds to expression ^ above if we choose U = °° (resp. L = 0). Formulae (1.1) and (1.2) in lH 
are used to define the corridor-realized variance in a way which treats the entrance of S into the corridor 
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differently from its exit from the corridor. This asymmetry is then exploited to obtain an approximate 
hedging strategy for linear payoffs on the corridor-realized variance. In this paper we opt for a symmetric 
treatment of the entrance and exit of S into and from the corridor [L, ?7], because this is in some sense more 
natural. It is however important to note that all the theorems and the algorithm proposed here do NOT 
depend in any significant way on this choice of definition. In other words for any reasonable modification of 
the definition in ^ (e.g. the one in the algorithm described in this section would still work. Note also 
that our algorithm will yield an approximate distribution of random variable ^ in the model S and therefore 
allows us to price any non-hnear payoff that depends on the corridor-realized variance. 

In the case the corridor-realized variance is monitored continuously (see |61), we can express it using 
quadratic variation as follows. Note first that since the map s max{L,mm{s,U}} can be expressed as 
a difference of two convex functions, Theorem 66 in [J 8 J implies that the process S = {St)t>o is again a 
semimartingale and therefore the corridor-realized variance 2^'^ (5), defined as the limit of the expression 
in (O as « tends to infinity, exists and equals 

(4) Qr^'^iS) = [log(S)]T-(logj) £ [l(o,L)(5,_)l(f;,.)(5,) + l(o,L)(5r)l(j/H(5r-)] 

V ^/ o<f<r 

by Theorem 4.47a in lfT3]| . Since we are assuming that the process S is cadlag the limit 5/_ '. — lim^ Sg exists 
almost surely for all t >0. The sum in (HJl, which corresponds to jumps of the asset price 5 over the corridor 
[L,?7], is almost surely finite by Theorem 4.47c in |[T3l . Note also that if L = (resp. U = °°) we find that 
Qy^ (S) (resp. Qy°°{S)) equals the quadratic variation of the semimartingale log(5) = (log(5'f))(>o because 
the process S cannot in these cases jump over the corridor. Our main task it to find an approximate law of 
the random variable Qj^ (S) which will allow us to price any derivative on the corridor-realized variance 
with terminal value 0(2^'^ (5)), where (p is a possibly non-hnear function. 

2.2. Markov chain X and its corridor-realized variance. Let us start by assuming that we are given a 
generator matrix ^ of a continuous-time Markov chain X = {Xt)t>o which approximates the generator of 
the Markov process S. The state-space of the Markov chain X is the set £ := {xq, . . . ,xn^\} C M+ with 
N elements, such that < xj for any integers < / < 7 < A'^ — 1. In Sections|4]and|5]we discuss briefly 
how to construct such approximate generators for Markov processes that are widely used in finance (i.e. 
diffusion processes with jumps.) Throughout the paper we will use the notation ^{x,y) = e'^^Cy for the 
elemetns of the matrix ^ , where x^y G E, vectors Cx^ey denote the corresponding standard basis vectors of 
and ' is transposition. 

The quantities of interest are the quadratic variation [log(X)] = ([log(X)],),>o and the corridor-realized 
vaiance = (2f '^(X)),>o processes which are for any maturity T defined by 

(5) [log(X)]r := lim £ (log 

f,"en„,,>i V 

(6) Qf{X) := [log(X)]r- (logy) £ [l(o,L)(X,_)l(j;,^)(X,) + l(o,L)(X,)l(j/,.)(X,_)] , 

V / ()</■< r 

where partitions n„, « G N, of [0, T] are as in ([U, the process X = {Xt)t>o is defined analogously with dJ]) by 
Xf := ma.x{L,min{Xf,U}} and := lim^y^Xs for any t > 0. Note that if we choose L < min{;c : x G £} 
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and U > max{x : x € £}, then the random variables in and Q coincide. We can therefore without loss 
of generality only consider the corridor-realized variance Q^^ {X). 

Since the process X is a finite-state Markov chain, the jumps of X arrive with bounded intensity and it is 
therefore clear that the following must hold 



(7) 



X 



t+ht 



Xt 



:o{At) for all x£[L,U]r\E. 



An analogous equality holds if Xt is ountside of the corridor [L,U]. Recall also that by definition a function 
f{At) is of the order o{At) (usually denoted by f{At) = o{At)) if and only if limA/\o/(AO /At = 0. Equal- 
ity ([7]l implies that the j-th instantaneous conditional moment of the corridor-realized variance Q^'^ {X) is 
given by 



Mj{x) 



(8) 



lim — E 

Af^O At 



yeE 



2; 



log 



log 



U 



Xt 

2; 



1a„ 



where the set Ay ^ C is defined as Aj/j, := ([0,L) x (?7,oo)) u {{U ,oo) x [0,L)) and for any x G £ we have 
X := max{L,min{x, ?7}}. 



2.3. The extension {X,l). The basic idea of this paper is to extend the markov chain X to a continuous- 
time Markov chain {X,l) = {Xt,It)t>o where the dynamics of the process / approximates well the dynamics 
of the corridor-realized variance Q^'^(X). Conditional on the path of the chain X, the process / will be a 
compound Poisson process with jump-intensity that is a function of the current state of X. The generator 
of {X,I) will be chosen in such a way that the first ^ G N infinitesimal moments of / and Q^'^ (X) coincide. 
The approximating chain / will start at (as does the process Q^'^ (X)) and gradually jump up its uniform 
state-space {0, a, .., a2C}, where a is a small positive constant and C is some fixed integer. 

The main computational tool in this paper is the well-known spectral decomposition for partial-circulant 
matrices (see Appendices A.2-A.4 in [IJ for the definition and the properties of the spectrum), which will 
be applied to the generator of the Markov chain {X,I). The geometry of the state-space {0,a,..,a2C} is 
therefore of fundamental importance because it allows the generator of (X,/) to be expressed as a partial- 
circulant matrix. As mentioned in the introduction, the main difference between the approach in the present 
paper and the algorithm in 111 is that here we take advantage of the full strength of the partial-circulant 
form of the generator of {X,I). This allows us to define the process / as a compound Poisson process 
with state-dependent intensity rather than just a Poisson process (which was the case in yj), without adding 
computational complexity. As we shall see in Section[6]this enables us to approximate the entire distribution 
of the corridor-realized variance and hence obtian much more accurate numerical results. 

Assuming that the process / can jump at most « G N states up from its current position in an infinitesimal 
amount of time, the dynamics of / are uniquely determined by the state-dependent intensities 



(9) 



A,-:£'^R+, where /g{1,...,?i} 
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and E is the state-space of the chain X. The generator of /, conditional on the event Xt = x, can therefore for 
any c,<i G {0, 1, ..,2C} be expressed as 

{Xj{x) ifd = c + i mod(2C+l) for some 7e{l,...,«}; 

-L1=ih{x) if d = c- 
otherwise. 

The dimension of the matrix (x : • , • ) is 2C + 1 for all x G £ and the identity d = c + j mod (2C + 1 ) means 
that the numbers d and c + j represent the same element in the additive group Z2C+1 ■ A key observation here 
is that the entries ^'{x : c,d) in the conditional generator depend on c and d solely through the difference 
d — c and hence the afore mentioned group structure makes the conditional generator into a circulant matrix 
(see Appendix for the definition of circulant matrices). 

This algebraic structure of the conditional generator ^\x : •,•) translates into aperiodic boundary con- 
dition for the process /. This is very undesirable because the process Q^-^ {X) we are trying to approximate 
clearly does not exhibit such features. We must therefore choose C large enough so that even if the chain 
/ is allowed to jump n steps up at any time, the probability that it oversteps the boundary is negligible (i.e. 
below machine precision). We will see in Section|6]that in practice C « 100 and « 30 is sufficient to avoid 
the boundary. Since our aim is to match the first k instantaneous moments, it is necessary to take n larger 
or equal to k. In applications this does not pose additional restrictions because, as we shall see in Section |6l 
k = 3 produces the desired results for jump-diffusions and ^ = 2 is already enough for continuous processes. 

The conditional generators given in (ITOl) can be used to specify the generator of the Markov chain (X , /) 
on the state-space £ x {0, a, ... , alC} as follows 

(11) ^ix,c;y,d) := ^ix,y)5e,d + (x : c,d)5,^y, 

where x,y £ E, c,d ^ {0,1,... ,2C} and 5.,. denotes the Kronecker delta function. The matrix is of the 
size N{2C+ 1) and has partial-circulant form. In other words we can express ^ in terms of blocks where 
each block is a square matrix of the size 2C + 1 and the blocks that intersects the diagonal of are equal to 
a sum of a circulant matrix and a scalar multiple of the identity matrix. All other blocks are scalar matrices. 
For the precise definition of partial-circulant matrices see Appendix 1X1 

We can now compute, using (ITOl ) and (fTTI) . the j-th instantaneous conditional moment of the process / as 
follows 

2C 

= Y^{ad-acy^'{x:c,d) 

d=0 

(12) = ajf^da^ix) 

d=\ 

for any x G £ and all integers c G {0, 1 , . . . , 2C} that satisfy the inequality c <2C — n, where n was introduced 
in Q. This inequality implies that the process / cannot jump to or above a2C (i.e. it cannot complete a 
full circle) in a very short time interval Af. Note also that it is through this inequality only that identity ([T2l) 
depends on the current level ac of the process /. 

Our main goal is to approximate the process {X,Q^'^ {X)), where corridor-realized variance Q^'^ {X) is 
defined in by the continuous-time Markov chain {X,I) with generator given by (fTTT) . We now match the 



hm — E 

Af->n Ait 



f+Af ■ 



Xf — x,It 



ac 
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first k instantaneous conditional moments of processes Q^'^ {X) and / using identities ([8]) and (fT2)) : 

n 

(13) ^ d^h{x) = Mj{x) for any xeE and 7 = 1 , . . . , fc. 

d=i 

In other words we must choose the intensity functions A,- (see Q) and the parameter a so that the system (fT3] ) 
is satisfied. The necessary requirement for the solution is that Xi{x) > for all x ^ E and all / = 1, . . . 
These inequalities can place non-trivial restrictions on the solution space and will be analysed in more detail 
in Sections |4] and [5] 

Another simple yet important observation that follows from (fT3l) is that, in order to match the first k 
instantaneous conditional moments of the corridor-realized variance Q^'^ {X), the size of the support of the 
jump distribution of the of Poisson processes with state-dependant intensity (i.e. n) must be at least k. From 
now on we assume that n> k. 

The pricing of volatility derivatives is done using the following theorem which yields a closed-form 
formula for the semingroup of the Markov chain (X,/). 

Theorem 2.1. Let 'S he the generator matrix of the Markov process {X,I) given by (II II ). Then for any t >0, 
x,y G E and d G {0,. .. , 2C} the equality holds 

F{X,=y,It = ad\Xo=x) = exp{t'^){x,0;y,d) 

1 2C 

(14) = ^^^Y.^e'P^'cMt^j){x,y), 

where i = \/— T, the scalars pj and the complex matrices for j = 0,... , 2C, are given by 

(15) ^jix,y) := ^{x,y) + 5,^yf^{e-'P^'-l)Xiix), 

(=1 

_ 271 
2C+1-'' 

Theorem 12. H is the main computational tool used in this paper which allows us to find in a semi-analytic 
form the semigroup of the chain {X,I) (if C 100 and = 70, the matrix contains more than 10^ 
elements). For a straightforward implementation of the algorithm in Matlab see |[T4l . It is clear that The- 
orem 12.11 generalizes equation (6) in [ll and that this generalization involves exactly the same number of 
matrix operations as the algorithm in HI. The only additional computations are the sums in ([T5] ). 

The proof of Theorem 12. 1 1 relies on the partial-circulant structure of the matrix 5^ given in (fTTI) . The 
argument follows precisely the same lines as the one that proved Theorem 3.1 in [Ij and will therefore not 
be given here (see Appendix A.5 in |[T1 for more details). 

Since the dynamics of the process {X,I) are assumed to be under a risk-neutral measure, the current value 
of any payoff that depends on the corridor-realized variance at fixed maturity can easily be obtained from the 
formulae in Theorem 12. II Furthermore the same algorithm yields the risk sensitivities Delta and Gamma of 
any derivative on the corridor-realized variance, without adding computational complexity. This is because 
the output of our scheme is a vector of values of the derivative in question conditional on the process X 
starting at each of the elements in its state-space. We should also note that forward-starting derivatives on 
the corridor-realized variance can be dealt with using the same algorithm because conditioning on the state 
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of a Markov chain at a future time requires only a single additional matrix-vector multiplication. Explicit 
calculations are obvious and are omitted (see HI for more details). 

3. Convergence 

In Section |2] we defined the Markov chain {X,I^) via its generator (ITTI ) that in some sense approximates 
the process {X,Q^--^ {X)), where Q^-^ {X) is the corridor-realized variance of X defined in Here 
denotes the process / from Section |2] which satisfies the instantaneous conditional moment restrictions, 
given by ([T3] ). up to order k. 

Notice that it follows directly from definition ^ that the process {X ,Q^'-^ {X)) is adapted to the natural 
filtration generated by the chain X and that its components X and Q^'^ {X) can only jump simultaneously. 
On the other hand note that the form of the generator of the chain {X,l'^), given by (ITTI ). implies that the 
components X and cannot both jump at the same time. It is also clear that the process l'^ is not adapted to 
the natural filtration of X. In this section our goal is to prove that, in spite of these differences, for any fixed 
time T the sequence of random variables {Ij)ken converges in distribution to the random variable Q^^ {X). 
In fact we have the following theorem which states that, for any bounded European payoff, the price of the 
corresponding derivative on the corridor-realized variance in the approximate model {X,I^) converges to the 
price of the same derivative in {X,Q^-^ {X)) as the number k of matched instantaneous conditional moments 
tends to infinity. 

Theorem 3.1. Let X be a continuous-time Markov chain with generator ^ as given in Section^ For each 
^ G N define a real number 

(16) at ■= ^max|(log^)^ : x,3;G£\{0}|, 

assume that n in ^ equals k and that there exist fiinctions X,f \ E ^ , / G { 1 , . . . , that solve the 
system of equations (1131 ). Let the continuous -time Markov chain (X,/*^) be given by generator (1111 ) where 
the integers Ck in (1101 ). which determine the size of the state-space of the process I^, are chosen in such a way 
that lim^^^oo (XuCk = °°- Then for any fixed time T > the sequence of random variables {I^)keN converges 
weakly to Q^^ {X). In other words for any bounded continuous function / : M — > R we have 

limE[/(4)|Xo]=E[/(e^'^(X))|Xo]. 

Before proving Theorem 13. II we note that the assumption on the existence of non-negative solutions of 
the system in ([T3] ) is not stringent and can be satisfies for any chain X by allowing « in (|9l) to take values 
larger than k. The restriction n = kin Theorem l3.1l is used because it simplifies the notation. 

Proof. Throughout this proof we will use the notation £, := Q^'^ {X) for any t G M+. By the Levy continuity 
theorem it is enough to prove that the equality holds 



hm E[exp(iM4)] = E[exp(iMr7')] for each m G M. 
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Let A? > be a small positive number and note that, by conditioning on the a-algebra generated by the 
process X up to and including time T — At and using the Markov property, we obtain the following repre- 
sentation 



(17) 



E[exp(iMrr)] = E[exp{±uLT-At)E[exjp{iu{'LT -'LT-At))\XT-At] 



E 



\j=o J- 



Aty\XT-At\ 



j=k+\ 



7! 



Af 



E 



7=1 



J! 



+ £ i^E[(rr-i:r-A.)^lXr-A.] 

j=k+l J- 



+ o{At) 



where Mj is defined in (H). By applying Markov property of {X,l'^), identity (fT2l) and condition ([T3l) . which 
holds by assumption for all j G {1, . . . ,k}, we obtain 



(18) 



E[exp(iM7; 



E 



(iuV 



e^-'r-.li+Atl^S^MjiXr^At) 



j=k+i 



T ~^T-At 



Xt- 



At,lT-At 



+ oiAt). 



It follows from ([8]) that there exists a positive constant G such that max{My(;c) : x E} <Gj for all j G N. 
Therefore we find that for a constant D := exp(MG) the following inequahty holds on the entire probability 
space 



(19) 



<D. 



Note also that D is independent of k and At. 

Definition ([T6l ) implies that /cttyt is a positive constant, say A, for each ^ G N. If we introduce a positive 
constant L := max{— (x,x) : x G £}, we obtain the following bound 



(20) 



EICLt -'LT-Aty\XT-At] < A^LAt + o{At) for each j G 



on the entire probability space. In order to find a similar bound for the process we first note that if follows 
from the linear equation ([T3] ) (for 7 = 1) and definition (fT6] ) that the inequalities 



k 

I 

d=\ 



Y^dX^{x)<kL for all keN,xeE 



must hold. Therefore ([12)) implies 



(21) 



E 



X 



T-At,iT-At 



< A'kLAt + o[At) for any j'GN 
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and any small time-step Af. We can now combine the estimates in ([TTl) . ([TSll . ([191) . (l20l) and (|2TI) to obtain 
the key bound 

(22) |E[exp(iMr7-)]-IE[exp(iM4)]| < \E[exp{iuLT-At)]-'K[exp{iul!}-_^)]\{\+AtD) 

+ L{k+\)At £ l^+o(AO. 

;■=<:+ 1 

The main idea of the proof of Theorem [TT] is to iterate the bound in (l22l ) times. This procedure yields 
the following estimates 

\E[exp{±uLT)]-E[exp{iu4)]\ < DAt{\ + AtoY^/^'^-^ + L{k + \)T £ + 

j=k+i J- 

Since the left-hand side of this inequality is independent of At, the inequality must hold in the limit as 
At \ 0. We therefore find 

(23) |E[exp(iMEr)]-E[exp(iM7^)]| <L(yt+l)r £ 

j=k+i j- 

The right-hand side of inequality (1231 ) clearly converges to zero as k tends to infinity. This concludes the 
proof of the theorem. |-j 

Theorem 13.11 implies that the prices of the volatility derivatives in the Markov chain model X can be 
approximated arbitrarily well using the method defined in Section [2l Our initial problem of approximating 
prices in the model based on a continuous-time Markov process S is by Theorem l3. 1 [ reduced to the question 
of the approximation of the law of S by the law of X. This can be achieved by a judicious choice for the 
generator matrix of the chain X. Since this is not the central topic of this paper we will not investigate the 
question further in this generality (see |(9l for numerous results on weak convergence of Markov processes). 
However in Sections |4] and |5] we are going to propose specific Markov chain approximations for diffusion 
and jump-diffusion processes respectively and study numerically the behaviour of the approximations for 
volatility derivatives in Section |6] 

4. The realized variance of a diffusion process 

Our task now is to apply the method described in Section |2] to approximate the dynamics of the corridor- 
realized variance of a diffusion processes. The first step is to approximate the diffusion process S which 
solves the stochastic differential equation (SDE) 

(24) ^ = ydt + a(^l-yw„ 

with measurable volatility function a : M+ —>■ M+, using a continuous-time Markov chain X. A possible way 
of achieving this is to use a generator for the chain X given by the following system of linear equations 

£if(x,y) = 0, 

yeE 

(25) Y^^{x,y){y-x) = yx, 

yeE 

2 



Yj£(x,y){y-xf = a x^ 
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for each x £ E.ln Appendix IB] we give an algorithm to define the state-space E of the chain X. In Section|6] 
we provide a numerical comparison for vanilla option prices in the CEV model, i.e. in the case a{s) := 
Oos^^^, and in the corresponding Markov chain model given by the approximation above. Note that a 
Markov chain approximation X of the diffusion 5 is in the spirit of HI and is by no means the only viable 
alternative. One could produce more accurate results by matching higher instantaneous moments of the two 
processes (see lfT6l for rates of convergence in some special cases). 

If the solution of SDE ((24l) is used as a model for the risky security under a risk-neutral measure we have 
to stipulate that 7 = r, where r is the prevailing risk-free rate in the economy. Therefore by the first two 
equations in system (|25] ) the vector in with cooridnates equal to the elements in the set E represents an 
eigenvector of the matrix for the eigenvalue 7. Hence we find 

(26) E [Xt \Xo = x]= 4 exp (t^) £ ye,, = e'^x, VxGE, 

where denotes the standard basis vector in that corresponds to the element x^E'm the natural ordering 
and the operation ' denotes transposition. Therefore, under the condition y = r, the market driven by the 
chain X will also have a correct risk-neutral drift. 

Once we define the chain X, the next task is to specify the process / that approximates well the corridor- 
realized variance 2^'^(X) defined in Q. As we shall see in Section |6l matching the first two moments 
(i.e. the case ^ = 2 in Section O is sufficient to approximate the corridor-realized variance dynamics of a 
diffusion processes. It is therefore necessary to take n>2, where n is the number of states the approximate 
variance process / can jump up by at any given time (see To have flexibility we use n much larger 
than 2, usually around 30. However in order to maintain the tractability of the solution of system ([T3] ) we 
make an additional assumption that the intensities A; in (|9l), for / = 2, . . . are all equal to a single intensity 
function A„ : £ ^ M+. To simplify the notation we introduce the symbol 

m 

(27) b"j''" := £ P , where j,n,m£N and m>n. 

l=n+l 

System ([T3] ) can in this case be solved explicitly as follows 

, , , aMi(x)bl'" -M2(x)b\'" 

(28) Ai(x) = ' , for any xeE, 



/^Q^ 1 ( \ M2{x)-aMx(x) 

(29) A„(x) = — -y— , for any xG£, 

a^(b^ —b{ ) 

where Mj{x) is given in ([8]). Since the functions are intensities, all the values they take must be 

non-negative. The formulae above imply that this is satisfied if and only if the following inequalities hold 

bl'" M2(x) 

(30) > > « for all xeE. 

b\'" - M,{x) - 

It is clear that the function x ^ M2{x) /M\{x), x ^ E, depends on the definition of the chain X both through 
the choice of the state-space E and the choice of the generator Figure [Tb]contains the plot of this function 
in the special case of the CEV model. Inequalities (l30l) are used to help us choose a feasible value for the 
parameter a which determines the geometry of the state-space of the process /. Note also that (l30l ) implies 
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that the larger the value of n is, the less restricted we are in choosing a. In Section [6] we will make these 
choices explicit for the CEV model. 

The generator of the approximate corridor-realized variance /, conditional on the chain X being at the 
level X, is in general given by Formula ([TOl l. In this particular case the non-zero matrix elements :c,d), 
c,d e {0,\,... , 2C}, are given by 

Xi{x) if d = {c + \) mod(2C+l); 

^'{x:c,d) := <^ X„{x) ifd = {c + i) mod (2C+ 1), / G {2, . ..,«}; 

— Ai(x) — — l)A„(x) if <i = c. 

This defines explicitly (via equations (|28] ) and (|29l )) the dynamics of the chain {X,I) if the original asset 
price process 5 is a diffusion. In Section |6] we will describe an implementation of this method when S 
follows a CEV process and study the behaviour of certain volatility derivatives in this model. 

5. The realized variance of a jump-diffusion 

In this section the task is to describe the algorithm for the pricing of volatility derivatives in jump-diffusion 
models. This will be achieved by an application of the algorithm from Section |2] with k = 3. In Section |6] 
we will investigate numerically the quality of this approximation. We start by describing a construction of 
the Markov chain which is used to approximate a jump-diffusion. 

5.1. Markov chain approximations for jump-diffusions. We will consider a class of processes with 
jumps that is obtained by subordination of diffusions. The prototype for such processes is the well-known 
variance gamma model defined in |15 |, which can be expressed as a time-changed Brownian motion with 
drift. 

A general way of building (possibly infinite-activity) jump-diffusion processes is by subordinating diffu- 
sions using a class of independent stochastic time changes. Such a time change is given by a non-decreasing 
stationary process (r,);>o with independent increments, which starts at zero, and is known as a subordinator. 
The law of {Tt)t>o is characterized by the Bernstein function 0(A), defined by the following identity 

(31) E[exp(-Ar,)] =exp(-0(A)O for any ?>0 and A G D, 

where D is an interval in M that contains the half-axis [0, °°). For example in the case of the variance gamma 
process, the Bernstein function is of the form 

(32) 0(A) = ^log(^l+A^). 

In this case (r,),>() is a gamma procesil] with characteristic function equal to IE[exp(iMrf)] = exp(— 0(— iw)?). 
Note that the set D in (|3TI) is in this case equal to (— /x/v,oo) (see |[T5l . equation (2)). This subordinator is 
used to construct the jump-diffusions in Section |6l 

Let 5 be a diffusion defined by the SDE in (l24l ). If we evaluate the process S at an independent subordi- 
nator (rf),>o, we obtain a Markov process with jumps (5'7;),>o. It was shown in fTTl that the semigroup of 
{ST,)t>o is generated by the unbounded differential operator := — 0(— ^), where $f denotes the generator 



^ The parameter /x is the mean rate, usually taken to be equal to one in order to ensure that E[r,] = / for all ? > 0, and V is the 
variance rate of (7;)/>o. 
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of the diffusion S. Similarly, if X is a continuous-time Markov chain with generator defined in the first 
paragraph of SectionlH the subordinated process {XT,)t>o is again a continuous-time Markov chain with the 
generator matrix := —<p{—^). We should stress here that it is possible to define rigorously the operator 
using the spectral decomposition of and the theorey of functional calculus (see 18], Chapter XIII, Sec- 
tion 5, Theorem 1). The matrix can be defined and calculated easily using the Jordan decomposition of 
the generator If the matrix ^ can be expressed in the diagonal form ^ = UAU^^, which is the case in 
any practical application (the set of matrices that cannot be diagonalised is of codimention one in the space 
of all matrices and therefore has Lebesgue measure zero), we can compute using the following formula 

(33) ^' = -U<p{-A)U-K 

Here <^'(— A) denotes a diagonal matrix with diagonal elements of the form 0(— A), where A runs over the 
spectrum of the generator 

Before using the described procedure to define the jump-diffusion process, we have to make sure that it 
has the correct drift under a risk-neutral measure. Recall that if the process S solves the SDE in (|24l) . then 
the identity E[Sf |So] = Soexp{tY) holds, where yis the drift parameter in (|24l) . Since the subordinator (r,),>o 
is independent of S, by conditioning on the random variable Tt, we find that under the pricing measure the 
following identity must hold 

5oexp(rt) =E[5r,|5o] =5oE[exp(7r,)] = 5oexp(-0(-7)O, 

where is the Bernstein function of the subordinator {Tt)t>o and r is the prevailing risk-free rate which is 
assumed to be constant. This will be satisfied if and only if r = — 0(— 7), which in case of the gamma sub- 
ordinator (i.e. when the function is given by (l32l) ) yields an explicit formula for the drift in equation (l24l) 

r ^ ^(i-exp(-r;)). 

Since formula (l26l ) holds for the chain X, tower property and the identity r = —(/)(— 7) imply 

E[Xt,\Xo] = XoE[exp (7r,)] =Xoexp{rt). 

Therefore the subordinated Markov chain {XT,)t>o can also be used as a model for a risky asset under the 
pricing measure. 

The construction of jump-diffusions described above is convenient because we can use the generator 
that was defined in Section IH and apply the Bernstein function from (l32l) to obtain the generator of the 
Markov chain that approximates the process (5'7;);>o. This accomplishes the first step in the approximation 
scheme outlined in the introduction. In Subsection 15.21 we develop an algorithm for computing the law of 
the relized variance of the approximating chain generated by It should be stressed that the algorithm 
in the next subsection does not depend on the procedure used to obtain the generator of the approximating 
chain. 

5.2. The algorithm. To simplify the notation let us assume that 5 is a jump-diffusion and that X is a 
coninuous-time Markov chain with generator ^ that is used to approximate the dynamics of the Markov 
process S. Since S has jumps it is no longer enough to use the algorithm from Section |2] with k = 2 (this 
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will become clear from the numerical results in Section [6]). In this subsection we give an account of how to 
apply our algorithm in the case k = 3. 

Assume we have chosen the spacing a and the constant C that uniquely determine the geometry of the 
state-space of the process / (see the paragraph preceding equation ^ for the definition of the state-space). 
Set the maximum jump size of / to be ma for some m G N. We now pick an integer n, such that 1 < « < m, 
and set the intensities that correspond to the jumps of the process / of sizes between 2a and na to equal 
Xn. Similarley we set the intensities for the jumps of sizes between {n + \)a and ma to be equal to A„,. 
This simplifying assumption makes it possible to describe the dynamics of / using only three functions 
A„, : E — > that give state-dependent intensities for jumping up by ia where / = 1, / G {2, . . . ,«}, 
/ e {n + 1, . . . ,m} respectively. In order to match k = 'i instantaneous conditional moments of the corridor- 
realized variance Q^^^ {X), these functions must by (fT3] ) satisfy the following system of equations 

Mj{x) 
ai ' 

the symbol Z?" '" is defined in (|27] ) and functions M,-, j = 1 , 2, 3, are given in dSjl. Gaussian elimination yields 
the explicit solution of the system 



/I 












1 








M2{x) 


WxGE, where Mj{x 


VI 


1 l.n 




\km{x)/ 


\m{x)) 





/Tl 1 n,m Yl i n,m\ i , \mi n,m i\,n,n.m\ /-jry ,n.m , n,m\ / 1 l,n i n,m ,\,n,n,m\ 

[Mjb^ —M\b^ )[b2 b{ —b{b2 j — yMjb^ — MiOj jip^^ b^ —b^'b^ ) 

{bY' - b''f){bYbT - bl-v^''") - {bi''" - bY"){bYbT - b\"b'r) ' 

(M3 - Ml ) (ft^''" - b^'"") - (Mj - Ml ) {b"/' - bl'"') 
^ ^ {m-M,){bY-bY) - {Ml -Mi){bl'"-b[") 

where all the identities are interpreted as functional equalites on the set E. It is clear from (|27] ) that the 
denominators in the above expressions satisfy the inequalities 

for suffciently large m (e.g. m > 10). This is because the term Zj"'™ dominates both expressions and has a 
negative coefficient in front of it. We therefore find that, if the functions Ai , A„ , A„, are to be positive, the 
following inequalities must be satisfied 

1 l.n, n.m , l,ni n.m , Im, umi il,n, n,m 

1 / N , . b^ b, — b, b^ , , o, —b, b^ 

(35) < a^Mi{x) + aM2{x)^ — \ — ^ M3 — \ — ^, 

^ ' '-y / ' ',n.mjl,n jl,njn,m -'^ 'jn.niii.n ;l,n;n.m' 

h h -bj' b^ bj b^ -b^' b^ 

yi.m yi,m jy","! j^n^m 

(36) > «^Mi (x) - aM2{x) J,,^ _ ^l,„ +M3(x) _ ^l,„ , 

(37) < a'M,{x)-aM2{x) \~\^ +M,{x) \~\^^ , 

^3' -^2 

for every x G These inequalities specify quadratic conditions on the spacing a (of the state-space of the 
process /) which have to be satisfied on the entire set E. 
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Note that inequality (1351) is always satisfied if the corresponding discriminant is negative. Alternatively 
if the discriminant is non-negative, then the real zeros of the corresponding parabola, denoted by cc{x),'a{x) 
and without loss of generality assumed to satisfy a{x) < a{x), exist and the conditions 

a < a{x) or a > a{x) Vx G £ 



must hold. Similar analysis can be applied to inequality (1371) . Inequality (1361) will always be violated if the 
discriminant is negative. This implies the following condition 



{b^ -b^ ){b^ -b{ ) M2{xY 



(38) ,,n,m jn,m-s(,n,m iri,m-s^ M^(^\2 Vx G 



which has to hold regardless of the choice of the spacing a. Even if condition (|38] | is satisfied we need to 
enforce the inequalities 

a(x)<a<a(x) Mx^E. 

In Table [U we summarise the conditions that need to hold for Ai (x) , A„ (x) , A„, (x) to be positive for any fixed 
element x G £. 





Discriminant > 


Restriction on a 


Ai(x) >0 


true 


a < a{x) or a > a(x) 




false 


none 


A„(x) > 


true 


a(x) < a < ■a(x) 




false 


A„ (x) cannot be positive 


A„, (x) > 


true 


a < a{x) or a > a(x) 




false 


none 



Table 1: Conditions on the discriminant and the real roots a(x),a(x), assumed to satisfy the relation a(x) < a(x), in this table 
refer to the parabolas that arise in inequalities i35\ . i36\ and ( I37t . These inequalities are equivalent to the conditions Xi (x) > 0, 
A„ (x) > and A,„ (x) > respectively. 



Having chosen the spacing a according to the conditions in Table \T\ we can use the formulae above to 
compute functions and A^. The conditional generator of the process /, defined in (flOl ). now takes the 

form 



^'{x:c,d) 



' Ai(x) ifc/ = (c + l) mod(2C+l); 
A„(x) if <i=(c + 5) mod (2C+ 1), 5 G {2, ...,«}; 
Am(x) ifJ=(c + 5) mod (2C+ 1), 5 G {«+ l,...,m}; 



with diagonal elements given by —X\ (x) — (« — 1 ) A„ (x) — {ni — n)Xm (x) and all other entries equal to zero. 
In Section[6]we are going to implement the algorithm described here for the variance gamma model and the 
subordinated CEV process. 
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6. Numerical results 

In this section we will perform a numerical study of the approximations given in the Sections |4] and |5l 
Subsection l6. ll gives an explicit construction of the approximating Markov chain X and compares the vanilla 
option prices with the ones in the original model S. Subsections 16.21 and [631 compare the algorithm for 
volatiUty derivatives described in this paper with a Monte Carlo simulation. 

6.1. Markov chain approximation. Let 5 be a Markov process that satisfies SDE (|24l) with the volatility 
function a : M+ M+ given by o{s) := Oqs^^^ and the drift / equal to the risk-free rate r (i.e. 5 is a CEV 
process). We generate the state-space E the algorithm in Appendix |B] and define find the generator matrix 
^ by solving the linear system in (l25l) . 

If the process 5 is a jump-diffusion of the form described in Subsection 15.11 (e.g. a variance gamma 
model or a CEV model subordinated by a gamma process), we obtain the generator for the chain X by 
applying formula (1331 ) to the generator defined in the previous paragraph, where the function <p (in (|33] )) 
is given by (l32l ). More preciselly if S is the subordinated CEV process, the drift 7 in (l24l) is given by the 
formula in (l34l ). If S is a variance gamma model we subordinate the geometric Brownian motion which 
solves SDE (l24l ) with the constant volatility function a{s) = Gq and the drift y = 6 + Gq /2 where 6 is the 
parameter in the vairance gamma model (see [15], equation (1)). The implementation in Matlab of this 
construction can be found in |[T4l . Note that the state-space of the Markov chain X in the diffusion and the 
jump-diffusion cases is of the same form (i.e. given by the algorithm in Appendix IB]). 

The numerical accuracy of these approximations is illustrated in Tables |2l [3] and |4] where the vanilla 
option prices in the Markov chain model X are compared with the prices in the original model S for the 
CEV process, the variance gamma model and the subordinated CEV model respectively. 





Markov chain X 


CEV: closed-form 


K\T 


0.5 


1 


2 


0.5 


1 


2 


80 


21.44% 


21.42% 


21.30% 


21.54% 


21.47% 


21.34% 


90 


20.55% 


20.57% 


20.46% 


20.68% 


20.62% 


20.49% 


100 


19.93% 


19.90% 


19.71% 


19.94% 


19.88% 


19.75% 


110 


19.37% 


19.19% 


19.11% 


19.28% 


19.22% 


19.10% 


120 


18.76% 


18.66% 


18.53% 


18.69% 


18.63% 


18.52% 



Table 2: Implied volatility in the CEV model. The maturity T varies from half a year to two years and the corresponding strikes 
are of the form Ke''^ , where K takes values between 80 and 120 and the risk-free rate equals r = 2%. The CEV process 5, with the 
current spot value 5o = 100, is given by ( 124b with the local volatility function CT equal to g{s) := CToi''^' and the drift y = r, 
where the volatility parameters are (Tq = 0.2, j8 = 0.3. The parameters for the non-uniform state-space of the chain X are N = 70 
and I = l,s = 100, J( = 700,^; = 50,^,, = 50 (see Appendix |B] for the definition of these parameters) and the generator of X is 
specified by system i25\ . The pricing in the Markov chain model is done using \39i and in the CEV model using a closed-form 
formula in 1121 . pages 562-563. The total computation time for all the option price in the table under the Markov chain model X is 
less than one tenth of a second on a standard PC with 1.6GHz Pentium-M processor and 1GB RAM. 

It is clear from Tables [2l [3] and |4] that the continuous-time Markov chain X approximates reasonably well 
the Markov process S on the level of European option prices. The pricing in the Markov chain model is 
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Markov chain X 


VG: FFT 


K\T 


0.5 


1 


2 


0.5 


1 


2 


80 


20.43% 


20.07% 


19.98% 


20.44% 


20.09% 


20.00% 


90 


19.91% 


19.89% 


19.93% 


19.95% 


19.94% 


19.96% 


100 


19.69% 


19.84% 


19.92% 


19.75% 


19.87% 


19.94% 


110 


19.85% 


19.89% 


19.93% 


19.82% 


19.88% 


19.93% 


120 


20.16% 


19.92% 


19.94% 


20.08% 


19.93% 


19.94% 



Table 3: Implied volatility in the variance gamma model. The strilces and maturities are as in Table|2] The process S, with the 
current spot value 5o = 100, is obtained by subordinating diffusion ( 124b with the constant volatility function f7(i') = CTo and the 
drift equal to y = + Oq /2, where 9 is given in [I5i , equation (1). The Bernstein function of the gamma subordinator is given 
in i32l . The risk-free rate is assumed to be r = 2%, the diffusion parameters take values cTq = 0.2, 6 = —0.04, and the jump 
parameters in i32i equal /x = 1 , v = 0.05. The parameters for the state-space of the chain X are N = 70 and / = i,s = 100, ii = 700, 
gl = 30, gu = 30 (see Appendix iBlfor the definition of these parameters). The total computation time for all the option price in the 
table under the Markov chain model is less than one tenth of a second. The Fourier inversion is performed using the algorithm 
in ||5J and takes approximately the same amount of time. All computations are performed on the same hardware as in Table|2] 





Markov chain X 


CEV with jumps: MC 


K\T 


0.5 


1 


2 


0.5 


1 


2 


80 


20.82% 


20.57% 


20.49% 


20.92% 


20.66% 


20.41% 


90 


20.08% 


20.10% 


20.11% 


20.16% 


20.12% 


20.19% 


100 


19.74% 


19.83% 


19.82% 


19.75% 


19.81% 


19.78% 


110 


19.66% 


19.61% 


19.56% 


19.64% 


19.58% 


19.53% 


120 


19.72% 


19.48% 


19.32% 


19.75% 


19.37% 


19.39% 



Table 4: Implied volatility in the CEV model subordinated by a gamma process. The strikes and maturities are as in Table|2] The 
process 5, with the current spot value So = 100, is obtained by subordinating diffusion l l24b with the volatility function 
a{s) = a^sP-^ (where cTq = 0.2, /3 = 0.7) and the drift given by i34\ (where the risk-free rate is r = 2% and the jump-parameters 
in l l32b equal fl = l.v = 0.05). The parameters for the state-space of the chain X are as in Table[3] The total computation time for 
all the option price in the table under the Markov chain model is less than one tenth of a second. The prices in the model S were 
computed using a Monte Carlo algorithm that first generates the paths of the gamma process {T,),^q (using the algorithm in fUl, 
page 144) and then, via an Euler scheme, generates paths of the process S. For the T = 2 years maturity, 10^ paths were generated 
in 200 seconds. All computations are performed on the same hardware as in Table|2] 

done by matrix exponentiation. The transition semigroup of the chain X is of the form 
(39) P(X, = y\XQ =x)= 4 &x^{t^)ey, where x,yeE, 

Cx^Cy are the corresponding vectors of the standard basis of and ' denotes transposition. For more 
details on this pricing algorithm see 0. The implied volatilities in the CEV, the variance gamma and 
the subordinated CEV model were obtained by a closed-form formula, a fast Fourier transform inversion 
algorithm and a Monte Carlo algorithm respectively. As mentioned earlier the quality of this approximation 
can be improved considerably, without increasing the size of the set E, by matching more than the first two 
instantaneous moments of the process S. 
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6.2. Volatility derivatives - the continuous case. The next task is to construct the process / defined in 
Section|2l obtain its law at a maturity T using Theorem 12. H and compaie it to the law of the random variable 
[log(5)]7- defined in ([Hi by pricing non-linear contracts. 

Let 5 be the CEV processs with the parameter values as in the caption of Tableland let X be the corre- 
sponding Markov chain, which is also described uniquely in the same caption. As described in Section |4] in 
tfiis case we use k = 2 (i.e. the process / matches the first and the second instantaneous conditional moments 
of the process [log(X)]7- defined in and hence define the state-dependent intensities in the conditional 
generator of by (l28l) and (l29l) . We still need to determine the values of the spacing a, the size (2C + 1) 
of the state-space of / and the largest possible jump-size an of the process / at any given time. 

The necessary and sufficien condition on parameters a and n is given by (l30l ). Figure [lb] contains the 
graph of the ratio in question x ^ M2{x) /M\{x), x ^ E, for the CEV model. The minimum of the ratio is 
0.000563, which can be used to define the value of a. The largest value of the ratio is approximately 0.019 
and hence « = 50 satisfies the first inequality in ([30l) . 

An important observation here is that Figure [Tbl only displays the values of the ratio M2{x)/M\{x) for 
xmEr\ [20,250]. The choices of a and n made above therefore satisfy condition ([30l ) only in this range 
(recall that in this case we have xq = 1 and xsg = 700). However this apparent violation of the condition 
in (l30l ) plays no role because the probability for the underlying process X to get below 20 or above 250 in 
2 years time is less than 10^^ (see Figure [Ta]). This intuitive statement is supproted by the quality of the 
approximation of the empirical distribution of [log(5)]7- by the distribution of Ij (see Figure[Tc]and Tabled. 

We now need to choose the size (2C + 1 ) of the state-space for the process /. The integer C is determined 
by the longest maturity that we are interested in, which in our case is 2 years. This is because we are using 
Theorem 12. II to find the joint law of the random variable {Xj,It) and must make sure that the process / 
does not complete the full circle during the time interval of length T (recall that the pricing algorithm based 
on Theorem 12.11 makes the assumption that the process / is on a circle). In other words we have to choose 
C so that the chain X accumulates much less than 2Ca of realized variance. In the example considered 
here it is sufficient to take C = 220, which makes the state-space {0, a, ... , a2C}, defined in the paragraph 
following ([S]), a uniform lattice in the interval between and 440 • 0.00056 = 0.246. Since the spacing a 
does not change with maturity, all that is needed to obtain the joint probability distribution of {Xt,It) for all 
maturities T £ {0.5, 1,2} is to diagonalize numerically the complex matrices ^j, j = 0, ...,2C, in ([T4l ) only 
once. The distribution of Ij, obtained as a marginal of the random vector {Xt,It), is plotted in Figure [Tel 
Note that the computational time required to obtain the law of It is therefore independent of maturity T . 

We now perform a numerical comparisons between our method for pricing volatility derivatives and a 
pricing algorithm based on a Monte Carlo simulation of the CEV model S. We generate 10^ paths of the 
process S using an Euler scheme and compute the empirical probability distribution of the realized variance 
[log(5')]7' based on that sample (see Figure [Tell. We also compute the variance swap, the volatility swap 
and the call option prices E [(Ir/r - {eKQf)+] , for G {80%, 100%, 120%}, where Kq := y^Epr/r] 
and Zj- denotes either Ij or [log(5')]7-. The prices and the computation times are documented in Table [5] A 
cursory inspection of the prices of non-linear payoff functions reveals that the method for ^ = 2 outperforms 
the algorithm proposed in HI, which corresponds to ^ = 1, without adding computational complexity since 
both algorithms require finding the spectrum of (2C + 1) complex matrices in ([T5] ). We will soon see that 
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Table 5: The prices of volatility derivatives in the CEV model S. The parameter values for the process 5 and the chain X are given 
in the caption of Table|2] The parameters for the process / are a = 0.00056, C = 220 for e { 1 ; 2} and n = 50 when k = 2 (recall 
from SectionUthat the parameter n controls the jumps of / strictly larger than a, which are note present if A: = 1). The variable Z7 
denotes either h or [log(5)]r and the call option price is E [(Zj- /T - {eKof)+] , for G {80%, 100%, 120%}, A'o := \/E[Sr/r]. 
An Euler scheme with a time-increment of one day is used to generate 10^ paths of the CEV process 5 and the sum in {TJ is used 
to obtain the empirical distribution of [log(5)]7- (see Figure [Tell and to evaluate the contingent claims in this table. The numbers in 
brackets are the standard errors in the Monte Carlo simulation. The computational time for the pricing of volatility derivatives 
using our algorithm is independent of the maturity T. All computations are performed on the same hardware as in TableEl 

the discrepancy between the algorithm in Q and the one proposed in the current paper is amplified in the 
presence of jumps. Note also that all three methods (k = 1,2 and the Monte Carlo method) agree in the case 
of linear payoffs. 

6.3. Volatility derivatives - the discontinuous case. In this subsection we will study numerically the 
behaviour of the law of random variables [log(5')]r and Qj^ (S), defined in ([T]) and (01) respectively, where S 
is a Markov process with jumps. Let 5 be a variance gamma or a subordinated CEV process with parameter 
values given in the captions of Tables [3] and |4] respectively. 

Since S has discontinuous trajectories we will have to match k = 3 instantaneous conditional moments 
when defining the process / in order to avoid large pricing errors for non-linear payoffs (see Tables |6] and [T] 
for the size of the errors when k = 2). We firts define the Markov chain X as described in Subsection l5. 1 l using 
the parameter values in the captions of Tables [3] and IH All computations in this subsection are performed 
using the implementation in |[T4l of our algorithm. 

Recall from Subsection 15.21 that in order to define the process / we need to set values for the integers 
I <n <m and the spacing a so that the intensity functions Ai,A„,Am are positive (Table [T] states explicit 
necessary and sufficient conditions for this to hold). Note that the inequality in (l38l) is necessary if A„ is 
to be positive. Figure l2c] contians the graph of the function x 1-^ 4Mi{x)Mi{x) / M2{x)'^ for x & E such that 
20 <x < 250 in the case of the variance gamma model. If we choose 



n:= 5 and m := 30, 



20 HARRY LO AND ALEKSANDAR MIJATOVIC 

then the left-hand side of the inequaUty in (l38l) equals 13.48, which is an upper bound for the ratio in 
Figure [2c] If S equals the subordinated CEV process, the graph of the function x h->- 4Mi{x)M^{x) /Mj^x)'^ 
takes a similar form and the same choice of ?i,m as above satisfies the inequality in (l38l) . 

The distance a between the consecutive points in the state-space of the process / has to be chosen so 
that the inequality a{x) < a < a{x) is satisfied for all x ^ E (see Table [B. Figure |2d] contains the graphs 
of the functions a, a over the state-space of X in the range 20 <x < 250 for the variance gamma model. 
The corresponding graphs in case of the subordianted CEV process are very similar and are not reported. It 
follows that by choosing 



a := 0.002 



we can ensure that all the conditions in the third row of Table [T] are met, both in the variance gamma and 
the subordinated CEV model, forx^E such that 20 <x < 250. It should be noted that it is impossible to 
find a single value of a that lies between the zeros a{x) and a{x) for all x ^ E for our specific choice of 
the chain X and its state-space. However not matching the instantaneous conditional moments of Ij and 
Q^'^ {X) outside of the interval [20,250] is in practice of little consequence because the probability that the 
chain X gets into this region (recall that the current spot level is assumed to be 100) before the maturity 
r = 2 is negligible (see Figure[2a]for the distribution of X in the case of variance gamma model). 

Once the parameters n,m and a have been determined, we use the explicit expressions for Ai, A„,A,„ on 
page[l4]to define the state dependent intensities of the process / for the states x G £ that satisfy 20 < x < 250. 
Outside of this region we choose the functions Xi,X„, A,„ : E —>■ to be constant. The choice of parameter 
C = 65 is, like in the previous subsection, determined by the longest maturity we are interested in (in our 
case this is T = 2). The laws of the realized variance [log(5')]r, for T G {0.5, 1,2}, in the variance gamma 
and the subordinated CEV model based on the approximation Ij are given in Figures l2bl and l3alrespectively. 
The prices of various payoffs on the realized variance in these two models are given in Tables |6] and |2l 

Observe that the time required to compute the distribution of Ij in the case of the continuous process S 
(see Table is larger than the time required to perform the equivalent task for the process with jumps (see 
Tables [6] and IT). From the point of view of the algorithm this difference arises because in the continuous 
case we have to use more points in the state-space of the process / since condition (l30l ) forces the choice 
of the smaller spacing a. In other words the quotient x i— > M2(x) /Mi (x) takes much smaller values if there 
are no jumps in the model S than if there are. It is intuitively clear from definition ^ that this ratio for the 
variance gamma (or the subordinated CEV) has a larger lower bound than the function in Figure[Tbl because 
in the the diffusion case the generator matrix is tridiagonal. 

Finally we apply our algorithm to computing the law of the corridor-realized variance Qj^ (S), where S 
is the subordinated CEV process and the corridor is given by L = 70 and U = 130. It is clear from Figure [3b] 
and the price of the square root payoff in Table [8] that the process / defined by matching k = 3 instantaneous 
moments of Q^'^ (S) approximates best the entire distribution of the corridor-realized variance. However, if 
one is interested only in the value of the corridor variance swap (i.e. a derivative with a payoff that is linear 
in Q^'^ (S)), Table [8] shows that it suffices to take k = I. 
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Table 6: The prices of volatility derivatives in the variance gamma model 5. The parameter values for the process 5 and the chain 
X are given in the caption of Table[3] The parameters for the process / are a = 0.002, C = 65 for k= 1,2,3. We choose n = 30 
when k = 2 and « = S.n? = 30 when = 3. The variable Ej- and the payoffs are as in Table|5] The algorithm in 1 1 1|, page 144, is 
used to generate 10^ paths of the VG process 5 and the sum in ([T} is used to obtain the empirical distribution of [log(5')]7- (see 
Figure [2bl ( and to evaluate the contingent claims in this table. The numbers in brackets are the standard errors in the Monte Carlo 
simulation. Note that the computational time for the pricing of volatility derivatives using the process / is independent of the 
maturity T . All computations are performed on the same hardware as in Table El ( see [14.1 for the source code in Matlab). 

7. Conclusion 

We proposed an algorithm for pricing and hedging volatility derivatives and derivatives on the corridor- 
realized variance in markets driven by Markov processes of dimension one. The scheme is based on an order 
k approximation of the corridor-realized variance process by a continuous-time Markov chain. We proved 
the weak convergence of our scheme as k tends to infinity and demonstrated with numerical examples that 
in practice it is sufficient to use ^ = 2 if the underlying Markov process is continuous and /: = 3 if the market 
model has jumps. 

There are two natural open questions related to this algorithm. First, it would be interesting to understand 
the precise rate of convergence in Theorem 13. ll both from the theoretical point of view and that of applica- 
tions. The second question is numerical in nature. As mentioned in the introduction, the algorithm described 
in this paper can be adapted to the case when the process 5 is a component of a two dimensional Markov 
process. The implementation of the algorithm in this case is hampered by the dimension of the generator of 
the approximating Markov chain, which would in this case be approximately 2000 (as opposed to 70, as in 
the examples of Section[6l). It would be interesting to understand the precise structure of this large generator 
matrix and perhaps exploit it to obtain an efficient algorithm for pricing volatility derivatives in the presence 
of stochastic volatility. 
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Table 7: The prices of volatility derivatives in the subordinated CEV model 5. The parameter values for the process S and the 
chain X are given in the caption of TablejH The parameters for the process /, the random variable Zj- and the payoffs of the 
volatility derivatives are as in Table[6] The algorithm described in the caption of Table|4]is used to generate 10^ paths of the 
process S and the sum in ([TJ is used to obtain the empirical distribution of [log(S)]7 (see Figure [Bat and to evaluate the contingent 
claims in this table. The numbers in brackets are the standard errors in the Monte Carlo simulation. Note that the computational 
time for the pricing of volatility derivatives using the process / is independent of the maturity T. All computations are performed 
on the same hardware as in Table|2](the code in 1 141 can easily be adapted to this model). 



CEV + jumps 


k moments 


Spectral: Ij 


MC: Q'f^iS) 


derivative \r 




0.5 


1 


2 


0.5 


1 


2 


corr-var swap 


1 


19.81% 


19.40% 


18.50% 


19.81% 


19.41% 


18.50% 


v/E[Er/r] 


2 


19.81% 


19.40% 


18.49% 


(0.051%) 


(0.050%) 


(0.048%) 




3 


19.81% 


19.40% 


18.50% 








corr-vol swap 


1 


19.59% 


19.22% 


18.25% 


19.12% 


19.03% 


18.19% 


Efv/Er/r] 


2 


19.18% 


18.97% 


18.08% 


(0.016%) 


(0.012%) 


(0.005%) 


3 


19.06% 


18.93% 


18.08% 








Time 




4s 


100s 


200s 


400s 



Table 8: Contingent claims on corridor-realized variance in the subordinated CEV model 5. The corridor is defined by L = 70 and 
U = 130. All parameter values are as in Table|7] The empirical distribution of iS) and the law of Ij for T G {0.5, 1,2} are 
given in Figure|3b] The Monte Carlo algorithm is as described in Table|4]and the numbers in brackets are the standard errors in the 
simulation. 
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Appendix A. Partial-circulant matrices 

A matrix C G M"^" is circulant if there exists a vector c G M" such that Cij = C(^i_j^ modn for iJ ^ 
{!,...,«}. The matrix C can always be diagonalised analytically, when viewed as a linear operator on the 
complex vector space C", as follows. For any r G {0, . . . ,« — 1} we have an eigenvalue and a corre- 
sponding eigenvector j^''^ (i.e. the equation Cy^''^ = Xrj^''^ holds for all r and the family of vectors y^''\ 
r G {0, — 1}, spans the whole of C") of the form 

A,= £c,e-'"'-*^ and = ^g-'fo for y g {o, 1}. 

It is interesting to note that the eigenvectors y'^^'\ r G {0, . . . ,?i — 1}, are independent of the circulant matrix 
C. For the proof of these statements see Appendix A in 111. 

Let A be a linear operator represented by a matrix in M*"^™ and let fi^*^' , for /: = 0, . . . , m — 1 , be a family 
of n-dimensional matrices with the following property: there exists an invertible matrix I] G C"^" such that 

U-^B^^^U = K^^^ , for all G {0, . . . ,m - 1}, 

where A^^) is a diagonal matrix in C"^". In other words this condition stipulates that the family of matrices 
B'^^^ can be simultaneously diagonalized by the transformation U . Therefore the columns of matrix U are 
eigenvectors of B^^"^ for all k between and m — \. 

Let us now define a large linear operator A, acting on a vector space of dimension mn, in the following 
way. Clearly the matrix A can be decomposed naturally into blocks of size « x «. Let A, y denote annxn 
matrix which represents the block in the i-th row and j'-th column of this decomposition. We now define the 
operator A as 

(40) An := B^') +A;;Im« and 

(41) Aij := Aijl^ii, for all /,7 G {1, . . . ,m} such that iy^j. 

The real numbers A,y are the entries of matrix A and I^n is the identity operator on E". We may now state 
our main definition. 

Definition. A matrix is termed partial-circulant if it admits a structural decomposition as in (l40l ) and (|4TI) 
for any matrix A G R'"^'" and a family of ^-dimensional circulant matrices B^'^^ for ^ = 0, . . . ,m — 1. 

For the spectral properties of partial-circulant matrices see Appendix A in 111. 

Appendix B. Non-uniform state-space of the Markov chain X 

The task here is to construct a non-uniform state-space for the Markov chain X, which was used in 
Section [6] to approximate the Markov process S. Recall that the state-space is a set of non-negative real 
numbers E = {xo,xi, . . . ,x/v_i} for some even integer G 2N. Recall that the elements of the set E, when 
viewed as a finite sequence, are strictly increasing. We first fix three real numbers l,s,u G M, such that 
I < s < u, that specify the boundaries of the lattice xq = I, x^^i = u and the stalling point of the chain 
x\N/2\ = s = So which coincides with the initial spot value in the model S. The function [•] : M ^ Z returns 
the smallest integer which is larger or equal than the argument. We next choose strictly positive parameter 
values gi,gu which control the granularity of the spacings between / and s and between s and u respectively. 



24 



HARRY LO AND ALEKSANDAR MIJATOVIC 



In Other words the larger gi (resp. gu) is, the more uniformly spaced the lattice is in the interval [l,s] (resp. 
[s,u\). The algorithm that constructs the lattice points is a slight modification of the algorithm in (19\ . page 
167, and can be described as follows. 

(1) Compute ci = arcsinh (^^) , C2 = arcsinh (^^^^ , Ni = \N/2] and = - (A^, + 1). 

(2) Define the lower part of the grid by the formula xt '■= s + g;sinh(ci (1 — k/Ni)) for k ^ {0, . . . ,Ni}. 
Note that xq = 1,xmi = s. 

(3) Define the upper part of the grid using the formula Xf^^j^k := s + ^„sinh(c2^/A'^«) for ^ G {0, . . . , A'^h}. 
Note that Xf^^i = u. 
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(a) The probability distribution function for the spot price Xj, 
with the maturity T equal to 0.5, 1 and 2 years, where X is the 
Markov chain used to approximate the CEV process 5. For a 
precise description of the process X see Subsection l6. II All 
relevant parameter values are given in the caption of Table|2] 
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(b) The function M2 (x) /Mi (x) , where x G £, in the CEV 
model. The minimum of this function, which equals 0.000563, 
determines the value of the spacing a by the second inequality 
in OOt . The maximum of the ratio, which is 0.019, determines 
the largest jump-size multiple n by the first inequality in ( I30t . 
All relevant parameter values for the CEV model and the 
accompanying chain X are given in the caption of Tabled 



PDF of the realized variance in tine CEV modei 
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X Monte Carlo [iog(S)]^(T=1/2) 
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(c) The empirical probability distribution of the realized variance [log(5)]7- of 
the CEV model 5, based on the Monte Carlo simulation described in 
Subsection l6.2l and the distribution of the random variable Ij, obtained from 
Theorem l2.1l for T e {0.5, 1,2}. For details on the definition of Ij see 
Sections[2]and|4] Note that the computational time required to obtain the law 
of It is independent of T (see caption of Table [3} . 



Figure 1 : CEV model 
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PDF of the spot price in the VG modei 



PDF of tile realized variance in the VG mode 




(a) The probability distribution function for the spot price 
Xj, with the maturity T equal to 0.5, 1 and 2 years, where X 
is the Markov chain used to approximate the variance 
gamma process. For a precise description of the process X 
see Subsection l6. II All relevant parameter values are given 
in the caption of Table[3] 
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(b) The empirical probability distribution of the realized variance 
[log(5)]7 in the VG model 5, based on the Monte Carlo 
simulation described in Subsection l6.3l and the distribution of the 
random variable Ij for T e {0.5, 1,2} matching k e {2,3} 
instantaneous moments. For details on Ij see Sections[2]and|5] 
Note that the computational time required to obtain the law of Ij 
is independent of T and that the quality of the approximation is 
greater for A: = 3 (see also Table|6]l. 



VG model - the ratio 4M {x)M (x)/M (x)"^ 



Solutions of the quadratic equation correspond to \ (x) in the VG model 



100 120 140 160 181 
X(x) 



200 220 240 



, for J e £ such that 



4Mi (x)M^(x) 

2Q<x< 250, in the variance gamma model. This function 



(c) The function x i 



appears in condition J38t of Subsection l5.2l The parameters 
of the chain X are given in the caption of Table[3] 
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(d) The functions x t-^ a(x) and x t-^ a(x), for x^E such that 
20 < X < 250, are the zeros of the quadratic in condition | |36I ) 
in the variance gamma model. As summarised in Table[T] in 
order to ensure that the intensity A„ [x) is positive, we must 
choose the value of the constant a to lie between the two 
curves for all x in the above range (see also Subsection l6.3t . 



Figure 2: Variance gamma model 
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PDF of the realized variance irn the subordinated CEV model 
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(a) The distribution of the realized variance in the subordinated CEV 
model. 

PDF of the corridor-realized variance in the subordinated CEV model 
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(b) The distribution of the corridor-realized variance in the 
subordinated CEV model. 



Figure 3: Figure[3al(resp.[3bt contains the empirical probability distribution of the realized variance [log(S)]7 (resp. 
corridor-realized variance Qj^ [S), where L = 70 and U = 130) in the subordinated CEV model 5, based on the Monte Carlo 
simulation described in the caption of Table|4](see also Subsection l6.3b . The distribution of the random variable Ij for the 
maturity T e {0.5, 1,2} matching k e {2,3} instantaneous moments is also plotted in both cases. For details on Ij see Sections|2] 
and|5] The computational time required to obtain the law of Ij is independent of T and the quality of the approximation improves 
drastically for k = 3 (see also Tables|7]and[8}. 



